/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
package technotheist.imgcog.net.predictor;

import technotheist.imgcog.base.Module;

/**
 *
 * @author Mitch
 */
public class PredictNet extends Module {

    protected void incrementT() {
        runner().T[0] += 1;
        runner().put(runner().T);
    }

    @Override
    protected PredictRunner runner() {
        return (PredictRunner) super.runner();
    }

    @Override
    protected void process() {
        throw new UnsupportedOperationException("Not supported yet."); //To change body of generated methods, choose Tools | Templates.
    }




    protected abstract class PredictRunner extends Runner {
        //

        protected final float[] outputErr;
        protected float[] inputErr;

        protected PredictRunner(PredictNetConfig config) {
            super(config);
            outputErr = new float[layerNSize];
            inputErr = new float[layer0Size];
            Q = config.getQLearnConstant();
            layerCnt = config.getLayerCount();
            X = new short[layerCnt];
            Y = new short[layerCnt];
            Z = new short[layerCnt];
            laySize = new int[layerCnt];
            layNode0 = new int[layerCnt];
            layN2P_Weight0 = new int[layerCnt];
            layP2L_Weight0 = new int[layerCnt];

            int nodeSum = 0;
            int n2pwSum = 0;
            int p2lwSum = 0;
            for (int lay = 0; lay < layerCnt; lay++) {
                X[lay] = config.getLayerXLen(lay);
                Y[lay] = config.getLayerYLen(lay);
                Z[lay] = config.getLayerZLen(lay);
                laySize[lay] = X[lay] * Y[lay] * Z[lay];
                layNode0[lay] = nodeSum;
                layN2P_Weight0[lay] = n2pwSum;
                layP2L_Weight0[lay] = p2lwSum;

                if (lay > 0) {
                    n2pwSum += laySize[lay] * laySize[lay - 1];
                }
                p2lwSum += laySize[lay] * LAT_CNT;
                nodeSum += laySize[lay];
            }
            //
            Pstore = new float[nodeSum];
            Lstore = new float[nodeSum];
            Mstore = new float[nodeSum];
            Fstore = new float[nodeSum];
            Nstore = new float[nodeSum];
            Perr = new float[nodeSum];
            Lerr = new float[nodeSum];
            Merr = new float[nodeSum];
            Ferr = new float[nodeSum];
            Nerr = new float[nodeSum];
            //
            //
            //
            N2P_Weight = new float[n2pwSum];
            for (int i = 0; i < n2pwSum; i++) {
                N2P_Weight[i] = config.getRandN2PWeight();
            }

            P2L_Weight = new float[p2lwSum];
            for (int i = 0; i < p2lwSum; i++) {
                P2L_Weight[i] = config.getRandP2LWeight();
            }
            //
            B2P_Weight = new float[nodeSum];
            //
            P2M_Weight = new float[nodeSum];
            P2F_Weight = new float[nodeSum];
            P2N_Weight = new float[nodeSum];
            //
            L2M_Weight = new float[nodeSum];
            L2F_Weight = new float[nodeSum];
            L2N_Weight = new float[nodeSum];
            //
            M2M_Weight = new float[nodeSum];
            M2F_Weight = new float[nodeSum];
            M2N_Weight = new float[nodeSum];
            //
            F2F_Weight = new float[nodeSum];
            F2N_Weight = new float[nodeSum];
            for (int i = 0; i < nodeSum; i++) {
                B2P_Weight[i] = config.getRandB2PWeight();
                P2M_Weight[i] = config.getRandP2MWeight();
                P2F_Weight[i] = config.getRandP2FWeight();
                P2N_Weight[i] = config.getRandP2NWeight();
                L2M_Weight[i] = config.getRandL2MWeight();
                L2F_Weight[i] = config.getRandL2FWeight();
                L2N_Weight[i] = config.getRandL2NWeight();
                M2M_Weight[i] = config.getRandM2MWeight();
                M2F_Weight[i] = config.getRandM2FWeight();
                M2N_Weight[i] = config.getRandM2NWeight();
                F2F_Weight[i] = config.getRandF2FWeight();
                F2N_Weight[i] = config.getRandF2NWeight();
            }

            //
            N2P_DW = new float[n2pwSum];
            P2L_DW = new float[p2lwSum];
            //
            B2P_DW = new float[nodeSum];
            //
            P2M_DW = new float[nodeSum];
            P2F_DW = new float[nodeSum];
            P2N_DW = new float[nodeSum];
            //
            L2M_DW = new float[nodeSum];
            L2F_DW = new float[nodeSum];
            L2N_DW = new float[nodeSum];
            //
            M2M_DW = new float[nodeSum];
            M2F_DW = new float[nodeSum];
            M2N_DW = new float[nodeSum];
            //
            F2F_DW = new float[nodeSum];
            F2N_DW = new float[nodeSum];
            //
            N2P_DWA = new float[n2pwSum];
            P2L_DWA = new float[p2lwSum];
            //
            B2P_DWA = new float[nodeSum];
            //
            P2M_DWA = new float[nodeSum];
            P2F_DWA = new float[nodeSum];
            P2N_DWA = new float[nodeSum];
            //
            L2M_DWA = new float[nodeSum];
            L2F_DWA = new float[nodeSum];
            L2N_DWA = new float[nodeSum];
            //
            M2M_DWA = new float[nodeSum];
            M2F_DWA = new float[nodeSum];
            M2N_DWA = new float[nodeSum];
            //
            F2F_DWA = new float[nodeSum];
            F2N_DWA = new float[nodeSum];
            //
            N2P_Mu = new float[n2pwSum];
            P2L_Mu = new float[p2lwSum];
            //
            B2P_Mu = new float[nodeSum];
            //
            P2M_Mu = new float[nodeSum];
            P2F_Mu = new float[nodeSum];
            P2N_Mu = new float[nodeSum];
            //
            L2M_Mu = new float[nodeSum];
            L2F_Mu = new float[nodeSum];
            L2N_Mu = new float[nodeSum];
            //
            M2M_Mu = new float[nodeSum];
            M2F_Mu = new float[nodeSum];
            M2N_Mu = new float[nodeSum];
            //
            F2F_Mu = new float[nodeSum];
            F2N_Mu = new float[nodeSum];
        }
        //
        //
        public final float Q;
        private final int[] T = new int[1];
        int layerCnt;
        private final short[] X;
        private final short[] Y;
        private final short[] Z;
        private final int[] laySize;
        private final int[] layNode0;
        //
        //**** **** **** ****
        //Values
        private final float[] Pstore;
        private final float[] Perr;
        private final float[] Lstore;
        private final float[] Lerr;
        private final float[] Mstore;
        private final float[] Merr;
        private final float[] Fstore;
        private final float[] Ferr;
        private final float[] Nstore;
        private final float[] Nerr;
        //
        //**** **** **** **** 
        //Connection weights
        private final int[] layN2P_Weight0;
        private final float[] N2P_Weight;
        private final float[] B2P_Weight;
        //
        private final int[] layP2L_Weight0;
        //
        private final float[] P2L_Weight;
        //
        private final float[] P2M_Weight;
        private final float[] L2M_Weight;
        private final float[] M2M_Weight;
        //
        private final float[] P2F_Weight;
        private final float[] L2F_Weight;
        private final float[] M2F_Weight;
        private final float[] F2F_Weight;
        //
        private final float[] P2N_Weight;
        private final float[] L2N_Weight;
        private final float[] M2N_Weight;
        private final float[] F2N_Weight;
        //---
        private final float[] N2P_DW;
        private final float[] B2P_DW;
        //
        private final float[] P2L_DW;
        //
        private final float[] P2M_DW;
        private final float[] L2M_DW;
        private final float[] M2M_DW;
        //
        private final float[] P2F_DW;
        private final float[] L2F_DW;
        private final float[] M2F_DW;
        private final float[] F2F_DW;
        //
        private final float[] P2N_DW;
        private final float[] L2N_DW;
        private final float[] M2N_DW;
        private final float[] F2N_DW;
        //---
        private final float[] N2P_DWA;
        private final float[] B2P_DWA;
        //
        private final float[] P2L_DWA;
        //
        private final float[] P2M_DWA;
        private final float[] L2M_DWA;
        private final float[] M2M_DWA;
        //
        private final float[] L2F_DWA;
        private final float[] P2F_DWA;
        private final float[] M2F_DWA;
        private final float[] F2F_DWA;
        //
        private final float[] L2N_DWA;
        private final float[] P2N_DWA;
        private final float[] M2N_DWA;
        private final float[] F2N_DWA;
        //---
        private final float[] N2P_Mu;
        private final float[] B2P_Mu;
        //
        private final float[] P2L_Mu;
        //
        private final float[] P2M_Mu;
        private final float[] L2M_Mu;
        private final float[] M2M_Mu;
        //
        private final float[] L2F_Mu;
        private final float[] P2F_Mu;
        private final float[] M2F_Mu;
        private final float[] F2F_Mu;
        //
        private final float[] L2N_Mu;
        private final float[] P2N_Mu;
        private final float[] M2N_Mu;
        private final float[] F2N_Mu;

        //
        void fwd(int layer, int node) {
            int pass = getPassId();
            if (pass == 0) {
                fwd_P(layer, node, layNode0[layer] + node);
            }
            if (pass == 1) {
                fwd_LMFN(layer, node, layNode0[layer] + node);
            }
        }

        void back(int layer, int node, int t) {
            int pass = getPassId();
            int ni = layNode0[layer] + node;
            float nX = Nstore[ni];
            float nE = 0;
            if (pass == 0) {
                nE = back_N(layer, node, ni, t);
            } else {
                nE = Nerr[ni];
            }
            float fX = Fstore[ni];
            float fE = Ferr[ni];
            if (pass == 0) {
                fE = back_F(layer, node, ni, t, nX, nE, fX, fE);
            }
            float mX = Mstore[ni];
            float mE = Merr[ni];
            if (pass == 0) {
                mE = back_M(layer, node, ni, t, nX, nE, fX, fE, mX, mE);
            }
            float lX = Lstore[ni];
            float lE = Lerr[ni];
            if (pass == 0) {
                lE = back_L(layer, node, ni, t, nX, nE, fX, fE, mX, mE);
            }
            if (pass == 1) {
                float pX = Pstore[ni];
                float pE = Perr[ni];
                pE = back_P(layer, node, ni, t, nX, nE, fX, fE, mX, mE);
            }
        }

        void fwd_P(int layer, int node, int index) {
            if (layer > 0) {
                float sum = 0;
                int i0 = layNode0[layer - 1];
                int w0 = layN2P_Weight0[layer] + laySize[layer - 1] * node;
                for (int i = 0; i < laySize[layer - 1]; i++) {
                    sum += Nstore[i0 + i] * N2P_Weight[w0 + i];
                }
                sum += B2P_Weight[index];
                Pstore[index] = biSigFunc(sum);
            } else {
                Pstore[index] = layer0Nodes[index];
            }
        }

        void fwd_LMFN(int layer, int node, int index) {
            fwd_L(layer, node);
            fwd_M(index);
            fwd_F(index);
            fwd_N(index);
        }

        void fwd_L(int layer, int node) {
            int nx = x3D(node, Y[layer], Z[2]);
            int ny = y3D(node, Y[layer], Z[2]);
            int nz = z3D(node, Z[2]);
            //
            int w = layP2L_Weight0[layer] + node * LAT_CNT;
            float sum = 0;
            for (int dx = -1; dx < 2; dx++) {
                for (int dy = -1; dy < 2; dy++) {
                    for (int dz = -1; dz < 2; dz++) {
                        sum += getP3D(layer, nx + dx, ny + dy, nz + dz)
                                * P2L_Weight[w];
                        w += 1;
                    }
                }
            }
            Lstore[layP2L_Weight0[layer] + node] = biSigFunc(sum);
        }

        void fwd_M(int index) {
            Mstore[index] =
                    Pstore[index] * P2M_Weight[index]
                    + Mstore[index] * M2M_Weight[index]
                    + Lstore[index] * L2M_Weight[index];
        }

        void fwd_F(int index) {
            Fstore[index] =
                    Pstore[index] * P2F_Weight[index]
                    + Mstore[index] * M2F_Weight[index]
                    + Lstore[index] * L2F_Weight[index]
                    + Fstore[index] * F2F_Weight[index];
        }

        void fwd_N(int index) {
            Nstore[index] = biSigFunc(
                    Pstore[index] * P2N_Weight[index]
                    + Mstore[index] * M2N_Weight[index]
                    + Lstore[index] * L2N_Weight[index]
                    + Fstore[index] * F2N_Weight[index]);
        }

        float back_N(int layer, int node, int ni, int t) {
            if (layer < layerCnt - 1) {
                float e = 0;
                int p0 = layNode0[layer + 1];
                int w0 = layN2P_Weight0[layer + 1];
                for (int p = 0; p < laySize[layer + 1]; p++) {
                    int wi = w0 + laySize[layer] * p + node;
                    e += backprop_N(p0 + p, ni, wi, t);
                }
                Nerr[ni] = e;
                return e;
            } else {
                float e = outputErr[node];
                Nerr[ni] = e;
                return e;
            }
        }

        float back_F(int layer, int node, int ni, int t, float nX, float nE, float fX, float fE) {
            float e = this.backprop_F(ni, t, nX, nE, fX, fE);
            Ferr[ni] = e;
            return e;
        }

        float back_M(int layer, int node, int ni, int t, float nX, float nE, float fX, float fE, float mX, float mE) {
            float e = this.backprop_M(ni, t, nX, nE, fX, fE, mX, mE);
            Merr[ni] = e;
            return e;
        }

        float back_L(int layer, int node, int ni, int t, float nX, float nE, float fX, float fE, float mX, float mE) {
            float e = this.backprop_L(ni, t, nX, nE, fX, fE, mX, mE);
            Lerr[ni] = e;
            return e;
        }

        float back_P(int layer, int node, int pi, int t, float nX, float nE, float fX, float fE, float mX, float mE) {
            float e = this.backprop_PInner(pi, t, nX, nE, fX, fE, mX, mE);
            int nx = x3D(pi, Y[layer], Z[layer]);
            int ny = y3D(pi, Y[layer], Z[layer]);
            int nz = z3D(pi, Z[layer]);
            int w0 = layP2L_Weight0[layer];
            int RX = X[layer];
            int RY = Y[layer];
            int RZ = Z[layer];
            for (int dx = -1; dx < 2; dx++) {
                for (int dy = -1; dy < 2; dy++) {
                    for (int dz = -1; dz < 2; dz++) {
                        //
                        int li = nodeIndex(layer, nx + dx, ny + dy, nz + dz);
                        int wi = w0 + (li - layNode0[layer]) * LAT_CNT + index3D(dx + 1, dy + 1, dz + 1, 3, 3);
                        e += backprop_PLateral(li, pi, wi, t);
                    }
                }
            }
            Perr[pi] = e;
            return e;
        }

        void back_B(int ni, int t, float pX, float pE) {
            backprop_B(ni, t, pX, pE);
        }

//        float pn_df(float fx) {
//            return 0.5f * (1f + fx) * (1f - fx);
//        }
        // P <-- [N_0 ... Nn]
        protected final float backprop_N(int pi, int ni, int wi, int t) {
            float x = 0;
            float e = 0;
            float delta = 0;
            //
            float dwa0 = 0;
            float mu0 = 0;
            float dwa1 = 0;
            float mu1 = 0;
            float dw = 0;
            //
            x = Pstore[pi];
            e = Perr[pi];
            delta = e * biSigDeriv(x);


            dwa0 = N2P_DWA[wi];
            mu0 = N2P_Mu[wi];
            //
            dw = mu0 * delta * x;
            dwa1 = nextDWA(dw, dwa0, t);
            mu1 = nextMu(mu0, dw, dwa0, dwa1);

            float ce = N2P_Weight[wi] * delta;
            N2P_Weight[wi] += dw;
            N2P_DWA[wi] = dwa1;
            N2P_Mu[wi] = mu1;
            return ce;
        }

        // N <-- F ; F <-- F
        protected final float backprop_F(int ni, int t, float nX, float nE, float fX, float fE) {
//            float x = 0;
//            float e = 0;
            float delta = 0;
            //
            float dwa0 = 0;
            float mu0 = 0;
            float dwa1 = 0;
            float mu1 = 0;
            float dw = 0;
            //
//            x = Nstore[ni];
//            e = Nerr[ni];
            delta = nE * biSigDeriv(nX);


            dwa0 = F2N_DWA[ni];
            mu0 = F2N_Mu[ni];
            //
            dw = mu0 * delta * nX;
            dwa1 = nextDWA(dw, dwa0, t);
            mu1 = nextMu(mu0, dw, dwa0, dwa1);

            float ce = F2N_Weight[ni] * delta;

            F2N_Weight[ni] += dw;
            F2N_DWA[ni] = dwa1;
            F2N_Mu[ni] = mu1;
            //
            //--- --- ---
//            x = Fstore[ni];
//            e = Ferr[ni];
            delta = fE;


            dwa0 = F2F_DWA[ni];
            mu0 = F2F_Mu[ni];
            //
            dw = mu0 * delta * fX;
            dwa1 = nextDWA(dw, dwa0, t);
            mu1 = nextMu(mu0, dw, dwa0, dwa1);

            ce += F2F_Weight[ni] * delta;

            F2F_Weight[ni] += dw;
            F2F_DWA[ni] = dwa1;
            F2F_Mu[ni] = mu1;
            return ce;
        }

        // F <-- M ; N <-- M ; M <-- M
        protected final float backprop_M(int ni, int t, float nX, float nE, float fX, float fE, float mX, float mE) {
//            float x = 0;
//            float e = 0;
            float delta = 0;
            //
            float dwa0 = 0;
            float mu0 = 0;
            float dwa1 = 0;
            float mu1 = 0;
            float dw = 0;
            //---- ----
            delta = nE * biSigDeriv(nX);

            dwa0 = M2N_DWA[ni];
            mu0 = M2N_Mu[ni];
            //
            dw = mu0 * delta * nX;
            dwa1 = nextDWA(dw, dwa0, t);
            mu1 = nextMu(mu0, dw, dwa0, dwa1);

            float ce = M2N_Weight[ni] * delta;

            M2N_Weight[ni] += dw;
            M2N_DWA[ni] = dwa1;
            M2N_Mu[ni] = mu1;
            //--- --- ---
            delta = mE;

            dwa0 = M2M_DWA[ni];
            mu0 = M2M_Mu[ni];
            //
            dw = mu0 * delta * mX;
            dwa1 = nextDWA(dw, dwa0, t);
            mu1 = nextMu(mu0, dw, dwa0, dwa1);

            ce += M2M_Weight[ni] * delta;

            M2M_Weight[ni] += dw;
            M2M_DWA[ni] = dwa1;
            M2M_Mu[ni] = mu1;
            //--- --- ---
            delta = fE;

            dwa0 = M2F_DWA[ni];
            mu0 = M2F_Mu[ni];
            //
            dw = mu0 * delta * fX;
            dwa1 = nextDWA(dw, dwa0, t);
            mu1 = nextMu(mu0, dw, dwa0, dwa1);

            ce += M2F_Weight[ni] * delta;

            M2F_Weight[ni] += dw;
            M2F_DWA[ni] = dwa1;
            M2F_Mu[ni] = mu1;
            return ce;
        }

        //M <-- P ; F <-- P ; N <-- P
        protected final float backprop_PInner(int ni, int t, float nX, float nE, float fX, float fE, float mX, float mE) {
            float delta = 0;
            //
            float dwa0 = 0;
            float mu0 = 0;
            float dwa1 = 0;
            float mu1 = 0;
            float dw = 0;
            //
            delta = nE * biSigDeriv(nX);

            dwa0 = P2N_DWA[ni];
            mu0 = P2N_Mu[ni];
            //
            dw = mu0 * delta * nX;
            dwa1 = nextDWA(dw, dwa0, t);
            mu1 = nextMu(mu0, dw, dwa0, dwa1);

            float ce = P2N_Weight[ni] * delta;

            P2N_Weight[ni] += dw;
            P2N_DWA[ni] = dwa1;
            P2N_Mu[ni] = mu1;
            //
            //--- --- ---
            delta = mE;


            dwa0 = P2M_DWA[ni];
            mu0 = P2M_Mu[ni];
            //
            dw = mu0 * delta * mX;
            dwa1 = nextDWA(dw, dwa0, t);
            mu1 = nextMu(mu0, dw, dwa0, dwa1);

            ce += P2M_Weight[ni] * delta;

            P2M_Weight[ni] += dw;
            P2M_DWA[ni] = dwa1;
            P2M_Mu[ni] = mu1;
            //
            //--- --- ---
            delta = fE;


            dwa0 = P2F_DWA[ni];
            mu0 = P2F_Mu[ni];
            //
            dw = mu0 * delta * fX;
            dwa1 = nextDWA(dw, dwa0, t);
            mu1 = nextMu(mu0, dw, dwa0, dwa1);

            ce += P2F_Weight[ni] * delta;

            P2F_Weight[ni] += dw;
            P2F_DWA[ni] = dwa1;
            P2F_Mu[ni] = mu1;
            return ce;
        }

        // L <-- [P_0 ... Pn]
        protected final float backprop_PLateral(int li, int pi, int wi, int t) {
            float x = 0;
            float e = 0;
            float delta = 0;
            //
            float dwa0 = 0;
            float mu0 = 0;
            float dwa1 = 0;
            float mu1 = 0;
            float dw = 0;
            //
            x = Lstore[li];
            e = Lerr[li];
            delta = e * biSigDeriv(x);


            dwa0 = P2L_DWA[wi];
            mu0 = P2L_Mu[wi];
            //
            dw = mu0 * delta * x;
            dwa1 = nextDWA(dw, dwa0, t);
            mu1 = nextMu(mu0, dw, dwa0, dwa1);

            float ce = P2L_Weight[wi] * delta;
            P2L_Weight[wi] += dw;
            P2L_DWA[wi] = dwa1;
            P2L_Mu[wi] = mu1;
            return ce;
        }
        // M <-- L ; N <-- L ; F <-- L 

        protected final float backprop_L(int ni, int t, float nX, float nE, float fX, float fE, float mX, float mE) {
//            float x = 0;
//            float e = 0;
            float delta = 0;
            //
            float dwa0 = 0;
            float mu0 = 0;
            float dwa1 = 0;
            float mu1 = 0;
            float dw = 0;
            //
            delta = nE * biSigFunc(nX);

            dwa0 = L2N_DWA[ni];
            mu0 = L2N_Mu[ni];
            //
            dw = mu0 * delta * nX;
            dwa1 = nextDWA(dw, dwa0, t);
            mu1 = nextMu(mu0, dw, dwa0, dwa1);

            float ce = L2N_Weight[ni] * delta;

            L2N_Weight[ni] += dw;
            L2N_DWA[ni] = dwa1;
            L2N_Mu[ni] = mu1;
            //
            //--- --- ---
            delta = mE;

            dwa0 = L2M_DWA[ni];
            mu0 = L2M_Mu[ni];
            //
            dw = mu0 * delta * mX;
            dwa1 = nextDWA(dw, dwa0, t);
            mu1 = nextMu(mu0, dw, dwa0, dwa1);

            ce += L2M_Weight[ni] * delta;

            L2M_Weight[ni] += dw;
            L2M_DWA[ni] = dwa1;
            L2M_Mu[ni] = mu1;
            //--- --- ---
            delta = fE;

            dwa0 = L2F_DWA[ni];
            mu0 = L2F_Mu[ni];
            //
            dw = mu0 * delta * fX;
            dwa1 = nextDWA(dw, dwa0, t);
            mu1 = nextMu(mu0, dw, dwa0, dwa1);

            ce += L2F_Weight[ni] * delta;

            L2F_Weight[ni] += dw;
            L2F_DWA[ni] = dwa1;
            L2F_Mu[ni] = mu1;
            return ce;
        }

        // M <-- L ; N <-- L ; F <-- L 
        protected final void backprop_B(int ni, int t, float pX, float pE) {
            float x = 0;
            float e = 0;
            float delta = 0;
            //
            float dwa0 = 0;
            float mu0 = 0;
            float dwa1 = 0;
            float mu1 = 0;
            float dw = 0;
            //
            delta = pE;

            dwa0 = B2P_DWA[ni];
            mu0 = B2P_Mu[ni];
            //
            dw = mu0 * delta * x;
            dwa1 = nextDWA(dw, dwa0, t);
            mu1 = nextMu(mu0, dw, dwa0, dwa1);

            B2P_Weight[ni] += dw;
            B2P_DWA[ni] = dwa1;
            B2P_Mu[ni] = mu1;
        }

        //*********************************************************
        protected float nextMu(float mu0, float dw, float dwa0, float dwa1) {
            return mu0 * max(0.5f, Q * dw * (dwa0 / (dwa1 * dwa1)));
        }

        protected float nextDWA(float dw, float dwa0, int t0) {
            float ex = 1f / (t0 + 1);
            return (dw * ex) + (dwa0 * (1f - ex));
        }
        //
        static final int LAT_CNT = 27;

        ///*******************************************************
        float getP3D(int layer, int x, int y, int z) {
            if (x >= 0 && x < X[layer] && y >= 0 && y < Y[layer] && z >= 0 && z < Z[layer]) {
                return Pstore[nodeIndex(layer, x, y, z)];
            } else {
                return Pstore[nodeIndex(layer,
                        (x + X[layer]) % X[layer],
                        (y + Y[layer]) % Y[layer],
                        (z + Z[layer]) % Z[layer])];
            }
        }

        int nodeIndex(int layer, int x, int y, int z) {
            return layNode0[layer]
                    + index3D(x, y, z, Y[layer], Z[layer]);
        }
    }
}